Collapse of differentially rotating supermassive stars: Post black hole formation 
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We investigate the collapse of differentially rotating supermassive stars (SMSs) by means of 3+1 
hydrodynamic simulations in general relativity. We particularly focus on the onset of collapse to 
understand the final outcome of collapsing SMSs. We find that the estimated ratio of the mass 
between the black hole (BH) and the surrounding disk from the equilibrium star is roughly the 
same as the results from numerical simulation. This suggests that the picture of axisymmetric 
collapse is adequate, in the absence of nonaxisymmetric instabilities, to illustrate the final state of 
the collapse. We also find that quasi-periodic gravitational waves continue to be emitted after the 
quasinormal mode frequency has decayed. We furthermore have found that when the newly formed 
BH is almost extreme Kerr, the amplitude of the quasi-periodic oscillation is enhanced during the 
late stages of the evolution. Geometrical features, shock waves, and instabilities of the fluid are 
suggested as a cause of this amplification behaviour. This alternative scenario for the collapse of 
differentially rotating SMSs might be observable by LISA. 

PACS numbers: 04.25. D-, 04.40. Dg, 97.10. Kc, 04.25. dg, 04.30.Db, 04.30. Tv 



I. INTRODUCTION 

There exists plenty of evidence that supermassive black 
holes (SMBHs) exist in the centre of galaxies, but their 
actual formation process has been a mystery for many 
decades Several different scenarios have been pro- 
posed, some based on stellar dynamics, others on gas 
hydrodynamics, and still others that combine these pro- 
cesses. There are three major routes to form a SMBH. 
One possibility is that the gas in the newly-forming 
galaxy does not break up into stellar-mass condensations, 
but instead forms a single superstar, and then undergoes 
complete gravitational collapse. Another possibility is 
that instead of forming a single superstar, black holes 
(BHs) might grow from smaller seeds, and then merge to 
form a SMBH. The final possibility is the runaway growth 
of a stellar mass BH to a SMBH by accretion. Here we 
consider the possibility of forming a SMBH from the col- 
lapse of a supermassive star (SMS). 

There are two categories of collapsing rotating SMSs 
based on their angular momentum distribution. One is 
the collapse of a uniformly rotating SMS. This happens 
when momentum transport is large, either through vis- 
cous turbulence or magnetic processes, which drives the 
star to rotate uniformly. In this case the path to the 
SMBH is the following. First the SMS contracts until the 
mass shedding limit is reached, conserving the angular 
momentum of the star. If the SMS contains sufficient an- 
gular momentum, the star evolves quasi-stationary along 
the mass-shedding limit, releasing the mass and angular 
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momentum. When the star reaches the post-Newtonian 
instability it begins to collapse. Baumgarte and Shapiro 
Q found that the star starts to collapse from the uni- 
versal configuration when the system contains sufficient 
angular momentum, independent of the history of the 
star. Saijo et al. [3] found that the collapse is coherent. 
That is, no significant bar instability occurs before BH 
formation in three dimensional post-Newtonian simula- 
tions. Shibata and Shapiro [4] found that a final Kerr BH 
of a/M w 0.75 contains 90% of the total rest-mass of the 
system, and that the disk around the BH contains 10% of 
that in an axisymmetric general relativistic simulation. 

The other category of collapsing rotating SMSs are dif- 
ferentially rotating. This happens when the viscous and 
the magnetic effects are small, which allows the star to ro- 
tate differentially. One of the representative scenarios for 
forming a differentially rotating star is as follows. First, 
a gas cloud gathers in an almost spherical configuration 
with some amount of angular momentum in the system. 
Next the almost spherical star contracts, conserving spe- 
cific angular momentum due to the lack of viscosity, to 
form a differentially rotating star, and possibly a disk at 
the end of the contraction J5j]. 

During the contraction of the differentially rotating 
SMS, prior to forming a supermassive disk, two possi- 
ble instabilities may arise that terminate the contrac- 
tion. One is the post-Newtonian gravitational instability, 
which leads the star to collapse dynamically. The other 
is the dynamical bar mode instability, which changes the 
angular momentum distribution of the star to form a bar, 
and possibly leads to the central core of the star collaps- 
ing to a BH due to the angular momentum loss. 

One of the primary observational missions for space- 
based detection of gravitational waves is the investigation 
of supermassive objects [6]. Since the Laser Interferome- 
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FIG. 1: Maximum of the rest mass density as a function 
of time. Solid, dashed, dotted, and dash-dotted lines denote 
models I, II, III, IV, respectively. 

ter Space Antenna (LISA) will have long arms (10 6 km), 
the detector will be most sensitive in the low frequency 
band (10 -4 ~ 10 _1 Hz). Potential sources of high sig- 
nal to noise events in this frequency range include quasi- 
periodic waves arising from nonaxisymmetric bars in col- 
lapsing SMSs and the inspiral of binary SMBHs (e.g. 
Ref. [7]). In addition, the nonspherical collapse of ro- 
tating SMSs to SMBHs could be a significant source of 
burst and quasi- normal gravitational waves (e.g. Ref. [21]). 
In this paper we track the collapse of a SMS by numerical 
simulation to investigate some of these possibilities. 

Here we focus on the post-Newtonian gravitational in- 
stability in differentially rotating SMSs. The aim of this 
work is to understand the final fate of the collapse of 
differentially rotating SMSs. There are two studies that 
address this issue. New and Shapiro [9[ studied the "zero- 
viscosity" equilibrium sequence in differentially rotating 
SMSs in Newtonian gravity. They showed that the star 
reaches the bifurcation point of dynamical bar instability 
before the disk forms by contraction. Saijo [10[ studied 
the gravitational collapse of differentially rotating SMSs 
in the conformally fiat approximation of general relativ- 
ity near the onset of radial instability. He showed that 
once the star collapses, the collapse is coherent and may 
form a SMBH with no substantial nonaxisymmetric de- 
formation due to instability. However the previous work 
may contain less angular momentum so that there might 
be no significant difference between the cases of uniform 
rotation and differential rotation. Therefore we now fo- 
cus on the case where the angular momentum is more 
contained in the star than the previous case, which po- 
tentially leads to rotational instabilities if they occur. 
In particular, we plan to answer the following questions. 
Does the BH form coherently? What are the features 
of the dynamics? Does the newly formed disk lead to 
instabilities? Can this system act as an efficient source 
of gravitational waves (GWs)? In order to answer these 



questions, three dimensional general relativist ic hydro- 
dynamics are desirable. 

The content of this paper is as follows. In Sec. HI 
we briefly explain the general relativistic hydrodynam- 
ics, especially the numerical tools we used to simulate 
the violent phenomena. In Sec. we introduce our 
findings of dynamical BHs, disk formation, and gravi- 
tational waves. Section IIVI is devoted to the summary 
of this paper. Throughout this paper, we use the ge- 
ometrized units with G = c = 1 and adopt Cartesian 
coordinates (x, y, z) with the coordinate time t. Note 
that Greek index takes (t,x,y,z), while Latin one takes 
(x,y,z). 



II. BASIC EQUATIONS 

In this section, we briefly describe the three- 
dimensional relativistic hydrodynamics in full general rel- 
ativity. 



A. The gravitational field equations 

We define the spatial projection tensor h^ v = g^ v + 
n^n v where g^ v is the spacetime metric, n M = 
(1/a, —f3 % /a) the unit normal to a spatial hypersurface, 
and where a and (3 1 are the lapse and shift. Firstly, we 
introduce a conformal factor ip = to the spatial three 
metric 7^ as 

lij = e- A Hj, (2.1) 

and choose the determinant of the conformally related 
metric 7^ to be unity, <j) = (hi7)/12. Secondly, we intro- 
duce a trace- free part of the extrinsic curvature A i3 - , and 
its conformally related tensor Aij, whose indices will be 
raised and lowered by the conformally related metric 7^ , 
as in Ref. 

A l3 = e~^A l3 = (k^ - \lijK^j , (2.2) 

where is the extrinsic curvature. Finally we introduce 
the conformal connection functions T % [12[ 

f^7^, = -^, (2.3) 

as a variable to describe the Ricci tensor of the confor- 
mally related metric. This step, the expansion of the 
Ricci tensor in terms of the metric based variables, plays 
a crucial role in maintaining a long time stable evolution 
of Einstein's field equations [111 ]. 

Therefore we evolve the spacetime with the 17 space- 
time associated variables ((/>, K, 7^, A i3 , P) as 
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TABLE I: Four different rotating equilibrium SMSs for evolution 



Model 


spacetime 


Rp /R e a 


^maxfe 


M c 


T/W d 


J/M 2e 


M/R f 


m disk 9 


(a/MY ati)tl 


I 


GR* 


0.450 


1.56 x 10 -5 


4.88 


0.108 


0.99 


2.56 x 10 -2 






I 


CF J 


0.450 


1.56 x 10" 5 


4.88 


0.108 


0.99 


2.56 x 10" 2 


0.044 


0.98 


II 


GR 


0.425 


1.56 x 10" 5 


5.07 


0.118 


1.03 


2.63 x 10" 2 






II 


CF 


0.425 


1.56 x 10" 5 


5.07 


0.118 


1.03 


2.63 x 10" 2 




>1 


III 


GR 


0.400 


1.56 x 10" 5 


5.31 


0.131 


1.07 


2.78 x 10" 2 






IV 


GR 


0.375 


1.56 x 10" 5 


5.75 


0.156 


1.10 


3.47 x 10" 2 







a Ratio of the polar proper radius to the equatorial proper radius 
b Maximum rest mass density 
c Gravitational mass 

d Ratio of the rotational kinetic energy to the gravitational binding 
energy 

e J: Total angular momentum 
fR: Circumferential radius 

ARatio of the estimated rest mass of the disk to the rest mass of 
the equilibrium star 
^Estimated Kerr parameter of the final hole 
General theory of relativity 
•^Conformally flat spacetime 



dt 



dt 



C 3 )T 



\aK, 
6 



-j lJ DiDja + a 
-2aAij, 



A ij A t3 + -K* + -( PH + S) 



e-^l-DiDja + a(i2y - S;,)] TF + otiKA^ - 2A il A l j ) 

c 

\*1 



2A^a + 2a (f^ - "-f^K - ^ + 6A^, 



dxi V dx l dx m 



3 7 dx lP 1 ' 



(2.4) 
(2.5) 
(2.6) 
(2.7) 

(2.8) 



where Cp denotes the Lie derivative along the shift f-driver [l5[ 

Ph = T^ v n^n v Si = T^n^h" and TF the trace- free part 
of the tensor. This set of equations when used to solve 
Einstein's field equations numerically is usually called the 
BSSN formalism. 

/ 8 3 \ 3 

We use as a slicing condition the generalised hyperbolic | j ^ — - a B\ (2.10) 



K-driver 0, [l| \dt d& J 4 

dt P dxo ) \dt P da* J V ' 

Ft " a = ~ a2f ^ K ~ ^ (2 ' U) 



where f(a) = 2/a, K$ = K t= o in order to avoid singular- 
ities. As for the shift, we use the generalised hyperbolic 



where r] is a damping parameter. 
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FIG. 2: Hamiltonian constraint (solid line) and Momentum 
constraints (dashed lines) for models I and II throughout the 
evolution. Note th at we normalised HC with 167rp ma x and 
Md with 87rp ma x y/M/R, where M is the gravitational mass, 
R the circumferential radius of the equilibrium star, respec- 
tively. 



B. The matter equations 

For a perfect fluid, the energy momentum tensor takes 
the form 

= p (l + £ + u»u v + Pg^, (2.12) 

where p is the rest-mass density, s the specific internal 
energy, P the pressure, the four- velocity. 
We adopt a T-law equation of state in the form 

P=(T-l)ps, (2.13) 

where T is the adiabatic index which we set to 4/3 in this 
paper, representing a SMS (the pressure is dominated by 
radiation). 

In the absence of thermal dissipation, Eq. (|2.13p . to- 
gether with the first law of thermodynamics, implies a 
polytropic equation of state 

P = ^p 1+1/n , (2.14) 

where n = 1/(T — 1) is the polytropic index and k a 
constant. 

From V ^T^ v = together with the continuity equa- 
tion, the flux conservative form of the relativistic con- 
tinuity, the relativistic energy and the relativistic Euler 
equations are [l6[ 

where the state vector U, the flux vectors T l ', and the 
source vectors S l are 



U 

S l 



[D,Si,T] T , 



0,7^ 



a 



D \ v l — — ] , Sj [ v l — — ] + PS),t v 1 -^- 



a 



9»i-r vil gsj),ay±- 



( T"°-S- In a -rt, 



(2.16) 
(2.17) 

(2.18) 



Note that (p, Vi, e) are the physical variables of the above 
equations, and the conserved quantities D, Si, r are 



D = pW, (2.19) 

W = olu\ (2.20) 

Si = phW 2 Vi, (2.21) 

E = phW 2 - P, (2.22) 

r = E-D, (2.23) 



where Vj = Uj/W, h = 1+s+P/p is the specific enthalpy. 
In the Newtonian limit, the above three physical variables 
coincide to the rest mass density, the flux density of the 
rest mass, and the energy of a unit volume of fluid. 
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FIG. 4: Centre of mass throughout the evolution. Solid and 
dash line denotes models I and II, respectively. 



The dynamical horizon mass of the hole Mbh can be 
computed once we have the angular momentum of the 
hole as 



1 

2Rs 



(2.25) 



Note that R s is the area radius of the hole. We also 
monitor the dimensionless Kerr parameter (J/M 2 )bh 
throughout the evolution. 



FIG. 3: Rest mass, ADM mass and the total angular mo- 
mentum throughout the evolution. We do not plot the above 
three quantities after the apparent horizon formation. Solid, 
dashed, dotted and dash-dotted line denote model I, II, III 
and IV, respectively. Since we use a different computational 
domain for computing a equilibrium and evolution quantities, 
the relative amount rate of ~ 0.5% appears in the ADM mass. 

C. Dynamical horizon 

A dynamical horizon is defined as the spacelike 
marginally trapped tube which is composed of future- 
marginally trapped surfaces on which the expansion as- 
sociated with the ingoing null direction is strictly nega- 
tive [l9[. If it is the outermost trapped surface, which is 
expected to be the case with all those considered here, 
then it is an apparent horizon. An angular momentum 
of the horizon Jbh can be written (see e.g. section III.B 
of Ref. HI) as 

Jbh = -^/ K^R^ds, (2.24) 

where W 1 is the outward directed spacelike normal to the 
horizon within the spacelike slice, and cp a is a rotational 
vector field on the horizon. This can be interpreted as 
the Komar angular momentum when (p a is a rotational 
Killing vector on the horizon. The code and method used 
to compute numerically an approximate Killing vector, 
should it exist, is described in Ref. [2l|. 



D. Gravitational waveform 

In order to monitor gravitational waves, we use the 
Newman-Penrose-Weyl scalar ^4 

^4 = C^ Xa k^fh u k x fh (T , (2.26) 

where C^xo- is the Weyl tensor, k^ the ingoing null vec- 
tor, ra M and ra M are the orthogonal spatial- null vectors 
of the four complex null tetrad (/ M , fc M , m M , m M ). More 
precisely, a set of null vectors are defined by introducing 
a orthonormal tetrad as 

k» = -L (efo - efo) , (2.28) 
m " = 7f( e (2)+ ie (3)). ( 2 - 2 9) 
™* = 7f ( e (2)- ie (3))- (2-30) 

Note that we choose the vector as the unit normal 
to the spatial hypersurface, the unit radial vector in 
spherical coordinates, and (e^, e^) the unit vectors in 
the angular directions. The Weyl scalar ^4 represents 
the outgoing gravitational waves at infinity 

V4 = h+-ih x , (2.31) 



6 



where h+ and h x are the two polarization modes 
(transverse-traceless condition) of the perturbed metric 
from fiat spacetime in spherical coordinate and q repre- 
sents the time derivative of the quantity q. If we put the 
observer sufficiently far from the source, the Weyl scalar 
^4 roughly represents the outgoing gravitational waves, 
ignoring the radiation back scattered by the curvature. 
Therefore we monitor the Weyl scalar ^4 to understand 
key features of the gravitational waves emitted by this 
system. 



for the collapsing stars. As in previous works we use 
a reconstruction-evolution method to evolve the matter, 
with PPM reconstruction [ItJ of the matter variables and 
Marquina's formula [l8[ to compute the inter-cell fluxes. 
Third order Runge-Kutta evolution in time is also used 
in our computation. 



III. NUMERICAL RESULTS 



A. Equilibrium Stars 



E. Treatment of the spacetime singularity from the 
newly formed black hole 

When a horizon is found in the simulation some of 
the matter within the horizon is excised from the do- 
main using the technique of Ref. [22]. The domain ex- 
cised is a sphere with coordinate radius half the min- 
imum coordinate radius of the horizon at the time of 
excision. It is allowed to grow during evolution but not 
to shrink (in coordinate radius). Although the spacetime 
is not excised, we introduce a strong artificial dissipation 
term in Eqs. (|2.4p - (|2.11|h combined with the standard 
gauge conditions, to avoid problems with the singularity 
[23[. In practice, we add a dissipation term —edx z dfq in 
Eqs. (|2.4p - (|2.1ip . where q is the quantity evolved by the 
particular evolution equation. We specify the dissipation 
parameter e ~ 0.2 for our simulation. 

In addition to this dissipation term, we have found it 
necessary to stabilize our evolutions in certain circum- 
stances by reducing the timestep on individual refine- 
ment levels. This is particularly true for the finest levels 
near the horizon where low density matter is orbiting at 
extreme speed. 



F. Computational codes and techniques 

We perform 3+1 hydrodynamic simulations in general 
relativity using CACTUS [24[ (gravitational physics), 
CARPET I25fl (mesh refinement of space and time), 
WHISKY (see Ref. [27| and references therein for a 
recent description) (general relativistic hydrodynamics). 
We set the outermost boundary of the computational grid 
for all direction as x max = 111 - 131M, imposing plane 
symmetry across the z = plane, and use 5-11 refine- 
ment levels to maintain the resolution at the central core 
region. We set the i-th refinement level of the uniform 
grid as U = 2 7 ~ i x 10, An = 2 5 ~* (i = 1, • • • , 11 for col- 
lapsing objects), where Li is the outer boundary of the 
grid, Ari is the stepsize of the spatial grid, respectively. 
The finest central resolution for the collapsing SMS is 
Ax = (3.08 - 3.21) x 10" 3 M. The Courant condition 
is set as At = C cr tAx, where C cr t = 0.05, 0.1, 0.2 for 
the finest three refinement levels for collapsing stars and 
Ccrt = 0-4 for the rest of the refinement levels, damping 
parameter of the shift condition rj = 3.91 - 4.05 M _1 



We construct the four differentially rotating SMSs for 
evolution. We use the polytropic equation of state of 
Eq. (|2.14p with n = 3, which represents the radiation 
pressure dominant, SMS sequence, and the relativistic 
"j-constant" rotation law as 



u u. 



= A 2 (n c -n), 



(3.1) 



where A(= R e A) represents the degree of differential 
rotation and has the dimension of length, Q c the cen- 
tral angular velocity, and Q the angular velocity. We 
choose A = 0.3 here. In the Newtonian limit (u* — > 1, 
u<p — > w 2 Vt where vo is the cylindrical radius of the star), 
this rotation law becomes 



A 2 Q r 



■A 2 ' 



(3.2) 



We use the technique of Ref. [28] to construct the ro- 
tating equilibrium stars. Our four differentially rotating 
equilibrium stars that will be evolved using the full T-law 
equation of state of Eq. (|2.13p are summarised in Table H 



B. Disk Mass from Equilibrium Configuration 

Before computing the collapse of the star, we can es- 
timate the final rest mass ratio between BH and its sur- 
rounding materials with the technique of Ref. [4]. The 
idea is the following. Suppose a star collapses axisym- 
metrically throughout the system. In this case the rest 
mass m* and the angular momentum j* along each cylin- 
der 



ra* = / dv^^gpu* , 

'Vcyl 



Jv rv : 



are, in the absence of viscosity, conserved throughout the 
evolution. Note that V cy \ is the three dimensional spa- 
tial volume along the cylinder. After the BH has formed, 
the material is swallowed at least up to the radius of 
the innermost stable circular orbit (ISCO) of the newly 
formed BH. Ignoring material from outside the ISCO we 
can approximate the rest mass of the BH, the angular 
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FIG. 5: Contours of the rest mass density in the equatorial plane for model I. Snapshots are plotted at values of (£, p ma x) 
= (a) (0, 1.56 x 10~ 5 ), (b) (671M, 5.04 x 10~ 2 ), (c) (755M, 5.28 x 10~ 5 ) and (d) (839M, 1.68 x 10~ 5 ), with a cutoff density 
of />cut = 1.56 x 10 -13 , respectively. The contour lines denote rest mass densities p/p ma x = I0~^ 15 ~^ d {i — 1, • • • , 14), where 
d — (log /?max — log p cu t)/15. Note that the apparent horizon exists after t — 671M, and the coordinate radius of the apparent 
horizon in the equatorial plane are Th. rz — (c) 0.283M and (d) 0.417M, respectively. 



momentum of the BH and the rest mass of the disk, de- 
fined as the mass outside the apparent horizon of the 
newly formed BH, as 

m™ = / dv^~gpu\ (3.3) 

Jv isco 

j* BH = / dv^pu'hu^ (3.4) 
Jv isco 

mf sk = [ dv^-gpu\ (3.5) 

^VW-V isco 

where ^ sco is the three dimensional spatial volume in- 
side the cylindrical radius of ISCO, Kvhi the whole three 
dimensional spatial volume, and V^hi — ^sco is the three 
dimensional spatial volume outside the ISCO cylindrical 
radius in each hypersurface. 



In practise, we compute the specific angular momen- 
tum at ISCO jisco- Note that j = hu^ is a monotonic 
growing quantity as a function of cylindrical radius, and 
that jisco is a function of the cylindrical rest mass ra* and 
the approximate, dimensionless Kerr parameter j*/m1. 
If there is a maximum in ji SCO as a function of the cylin- 
drical rest mass (cylindrical radius), the quantities at the 
maximum (cylindrical rest mass and specific angular mo- 
mentum) correspond to the quantity of the newly formed 
BH if the equilibrium star collapses. After the determi- 
nation of the cylindrical rest mass of the BH, we can 
compute that of the disk from the total rest mass of the 
system. We can also compute the approximate Kerr pa- 
rameter from the quantities at the maximum. To sim- 
plify the calculation our estimates of the ratio between 
disk mass and the BH mass and the Kerr parameter are 
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FIG. 6: Same as Fig. [5] but in the meridional plane. 



computed in conformally flat spacetime (using the meth- 
ods of e.g. Ref. [10[) and are summarised in Table H 



C. Collapse of Differentially Rotating 
Supermassive Stars 

We first investigate the onset of collapse by evolving 
the four differentially rotating equilibrium stars of Ta- 
ble H We choose the z-axis as the rotational one of the 
equilibrium star. Since we use the polytropic equation of 
state P = KpQ (P: pressure, n: constant, po: rest mass 
density, T: adiabatic exponent) when constructing initial 
data sets, all physical quantities are rescalable in terms of 
n. Therefore, we represent all physical quantities nondi- 
mensionally in this paper, which is equivalent to setting 
n = 1. To trigger collapse, we deplete pressure by 1%. 
Checking the maximum rest mass density of the rotating 
stars throughout the evolution, we conclude that models 
I and II are radially unstable, while models III and IV 
are stable (Fig. [1]). 

Since we use a free evolution scheme for the 17 vari- 
ables associated with the spacetime, the constraint equa- 
tions of Einstein's field equations are not automatically 
satisfied at each timestep. Therefore we can monitor the 
Hamiltonian constraint and Momentum constraints as a 
check on the accuracy of the simulation. We define them 
as 



HC 



R 



K — 167T/9H, 



ViK - 8wSi, 



(3.6) 
(3.7) 



and show them for models I and II in Fig. [2j Since these 
quantities should be exactly zero sense, we normalise the 
Hamiltonian constraint by 167 rp max and the Momentum 
constraints by 8irp mdiX ^ M/ R, where p max is the maxi- 
mum rest mass density at each timestep. Both of the 
normalisation quantities represent the typical size of the 
source term of the constraints (maximum rest-mass den- 
sity and collapsing momentum density). Since all con- 
straint quantities are conserved within ~1% from Fig. [21 
the constraints seem to be satisfied at a reasonable level 
throughout the evolution. 

We also monitor the rest mass Mo, ADM mass M a d m , 
and the total angular momentum J 



M 



AL 



adn 



dvpu 1 

1 

~2tt 



■/ 

1 

8^ 

-/ 



[(phW 2 - P)e 5( f + —e^AijA^], 

l07T 



dve**[xS y -yS x + 



16tt 



Aij(xd y - yd x )j^ 



+ ^(xd v -yd x )K k k ], 

throughout the evolution in Fig. [3l We only monitor 
these three quantities until the time when the appar- 
ent horizon forms, since our excision method significantly 
complicates the computation of the required volume in- 
tegrals after horizon formation. From Fig. O Mo is con- 
served to within 0.5%, M a d m to within 1%, and J to 
within 2%. Therefore, we can conclude that all three 
quantities are satisfactory conserved for our case. 

We also monitor the centre of mass condition in Fig. [H 
In practise we monitor 



M x 



Jdv (hal ^px 



M, 



(half) 



M, 



(half) 



(3.8) 
(3.9) 



where Mq a is the rest mass ofz > 0, dv( haU ^ the spatial 
volume of z > 0, since we adopt planar symmetry across 
z = 0. Since our finest spatial gridsize is « 3 x 10 _3 M, 
the centre of mass motion is contained within a few finest 
grids cells from the origin. Therefore, we can conclude 
that centre of mass condition is satisfied for our case. 

We show contours of the rest mass density in Figs. [5] 
(equatorial plane) and [6] (meridional plane) for model I 
and in Figs. (equatorial plane) and [8] (meridional plane) 
for model II. As we can easily understand from these fig- 
ures, the radially unstable stars first collapse to form a 
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FIG. 7: Contours of the rest mass density in the equatorial plane for model II. Snapshots are plotted at values of (£, p m ax) 
= (a) (0, 1.56 x 10" 5 ), (b) (768M, 5.07 x 10" 2 ), (c) (849M, 1.05 x 10" 4 ) and (d) (930M, 6.12 x 10" 5 ), with a cutoff density 
of p cut = 1.56 x 10" 



respectively. The contour lines denote rest mass densities p/p n 



10 



-(15-i)d 



(i 



, 14), where 



d = (log /?max — log p C ut)/15. Note that the apparent horizon exists after t = 770M, and the coordinate radius of the apparent 
horizon in the equatorial plane are r nr z — (c) 0.170M and (d) 0.175M, respectively. 



BEL After that, there is a significant amount of mass ejec- 
tion due to the large amount of angular momentum of the 
hole. It is indeed like a flare. Finally the system settles 
down to a quasi-stationary state; central BH with a disk. 
Fluid elements in the equilibrium star are balanced by its 
self gravity, the pressure gradient and centrifugal force. 
When the BH forms, some fluid elements act as a parti- 
cle which is free from the interaction. In this case, the 
material can spread out to a larger radius than the equi- 
librium radius of the star. Note that there seem to have 
m = 4 structure in the snapshots of I (b) (c) (Fig. [5]) and 
II (c) (Fig. [7j). We believe that this is due to the reduced 
spatial computational resolution of the outer regions of 
the star (See Ref. [10] for a different resolution case). 

Next we trace the mass and the angular momentum of 
the newly formed BH throughout the evolution using the 
dynamical horizon techniques outlined in Sec. Ill CI We 



monitor the gravitational mass, total angular momentum 
and the Kerr parameter of the newly formed BH through- 
out the evolution as shown in Fig. [9j The BH mass, the 
spin and the Kerr parameter increase monotonically after 
the BH has formed, by swallowing much of the surround- 
ing material. This stage lasts roughly until all of the 
matter located inside the radius of the innermost stable 
circular orbit of the final BH is swallowed. 

We have also confirmed that the estimated mass and 
spin of the BH from the equilibrium configuration of the 
collapsing SMS are in good agreement with the results 
from the dynamics. Although we calculated the disk 
mass from the equilibrium configuration in conformally 
flat spacetime, there is a little difference in the physical 
quantities between full general relativistic spacetime and 
conformally flat spacetime (see Table [I]). For instance, 
the estimated Kerr parameter from the equilibrium star 
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FIG. 8: Same as Fig. [7] but in the meridional plane. 



of model I is 0.98 (Table HJ), while the result of numerical 
simulation is « 0.97 (Fig. [9]). Also the ratio between the 
estimated rest mass of the disk and the rest mass of the 
equilibrium star of model I is 0.044 (Table [J), while the 
result of numerical simulation is ~ 0.05 (Fig. [TO]). Note 
that we define the rest mass of the disk as the rest mass 
outside the apparent horizon in each timeslice. 

We furthermore study the formation of a massive disk 
from the collapse of differentially rotating SMSs. We 
trace the rest mass of the disk for models I and II 
(Fig. [TO]) . The rest mass of the disk monotonically de- 
creases once the BH has formed, since the newly formed 
BH grows monotonically by swallowing the surrounding 
materials. One noticeable feature in Fig. [10] is that there 
is a plateau at the final stage of model II. This indicates 
that there is a strong angular momentum barrier that 
a fluid fragments are prohibited to fall into a hole. We 
also illustrate the snapshot of the rest mass density in 
the meridional plane of model II (Fig. [8] II [d]). The max- 
imum of the rest mass density is located around r « 2M 
in coordinate units. 

We also monitor the diagnostics of the m = 1 and 
m = 2 Fourier modes throughout the collapse in Figs.fTTl 
and [T2l (for model I) and Figs. [T3l and [TTI (for model II). 
We define the diagnostics C m as 



°d<p, 



(3.10) 
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FIG. 9: Gravitational mass (Mbh), total angular momentum 
(«7bh) and Kerr parameter ((J/M 2 )bh) of a newly formed BH 
as a function of time. Solid and dashed line represent models 
I and II, respectively. 




with a normalisation of Co, a mean density of the ring. In 
order to understand the characteristic frequency, we also 
show the spectrum of the ring diagnostics |.F[p m ](cc;)| 2 
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FIG. 10: Disk mass as a function of time. We defined the 
disk mass as the rest mass outside the apparent horizon of 
the newly formed BH. Solid and dashed line represent models 
I and II, respectively. 



using the Fourier transformation as 

FlPm](u>) = ^ dt pe^-^dip, (3.11) 

where t^ n is the time where we terminate our simulation. 
We show the spectra of the m = 1 and m = 2 diagnostics 
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FIG. 11: m = 1 and m — 2 diagnostics of the ring in the equa- 
torial plane throughout the collapse for model I. We monitor 
the diagnostics at the coordinate radius r = 1.02M (upper 
two panels) and 2.05M (lower two panels), respectively. Note 
that the solid line and the dashed line denote m — 1 and 
m = 2. The coordinate equatorial radius of the horizon at 
t 854M is r = 0.464M. 

in Figs. [T5l and [T6l 

The diagnostics show that non-axisymmetries (that 
would be triggered at least at the level of finite- 
differencing error through the simulation) do not grow 
to meaningful levels near the centre of the system un- 
til well past horizon formation. When these non- 
axisymmetries appear, they grow rapidly before saturat- 
ing. From the growth of the coefficients at larger radius, 
shown in Figs. [12] and HH it seems clear that the non- 
axisymmetries are triggered in the interior of the disk 
and propagate outwards. 



D. Gravitational radiation 

Finally we investigate the gravitational waveform from 
the collapsing object. We introduce the Weyl scalar ^ 4 
to study the outgoing gravitational waves. If we put the 
observer sufficiently far from the source, the Weyl scalar 
^4 roughly represents the outgoing gravitational waves, 



FIG. 12: Same as Fig. [TT] but we monitor the diagnostics 
at the coordinate radius r = 4.1 M (upper two panels) and 
10. 2M (lower two panels), respectively. 



since the radiation back scattered by the curvature be- 
haves as r -3 , where r is the distance from the source. 
We observe the waveform (the real component of ^4) 
along the x-axis in the equatorial plane at three differ- 
ent coordinate locations with radius r « 50M, « 65M, 
« 100M for models I and II (Figs. 13 and [[8]). Note that 
the equatorial coordinate radius of the equilibrium star 
is r = 38. 1M and 37. 0M, and the outermost boundary 
of the spatial grid is 131M and 126M, for models I and 
II. 

We find that the waveform contains three different 
stages. The first stage is the burst. This occurs around 
the time when the apparent horizon of the SMS forms. 
The dominant contribution of the burst comes from the 
axisymmetric mode due to collapse. The second stage is 
the quasinormal ringing of the newly formed BH. Since 
the dominant frequency in the spectrum (Figs. [19] and 
[20]) is Muj « 0.4, the dominant contribution is the ax- 
isymmetric one, using the fact that quasinormal mode 
frequency of / = 2, m = is Muj « 0.40 - 0.42 for 
a BH with Kerr parameter a > 0.8M [29], where /, m 
denotes the indices of spin — 2-weighted spheroidal har- 
monics. The final stage has a quasi-stationary wave. The 
amplitude of the waveform for model II seems to increase 
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FIG. 13: m = 1 and m = 2 diagnostics throughout the 
collapse for models II. We monitor the diagnostics at r = 
0.197M (upper two panels) and 0.987M (lower two panels), 
respectively. Note that the solid line and the dashed line 
denote m — 1 and m = 2. The coordinate equatorial radius 
of the horizon at t = 940M is r = 0.177M. 



at late times. The feature corresponds to several peaks 
in the spectrum (Fig. [20]h especially in the regime of 
Moo < 0.5. For model I, the quasi-stationary waveform 
remains for at least At ~ 70M (Fig. [20]) . In this case 
the spectra does not contain the peaks in the regime 
of Muo < 0.5. When we reduce the integration time of 
the quasi-periodic waves (comparing the three panels in 
Fig. [20]) , the spectrum becomes much smoother in the 
frequency regime Muj > 0.5. Therefore the frequency 
region of Moo > 0.5 plays a role in the quasi-periodic 
waveform. Although we do not know the origin of the 
late-time quasi-periodic waves, we can discuss possible 
candidates and their likelihood. 

Firstly we distinguish between the different possibili- 
ties for the amplification of the waves based on the dom- 
inant type of behaviour. The most obvious possibility 
would be strong matter dynamics perturbing the essen- 
tially stationary Kerr spacetime given by the dominant 
central BH. The second possibility is a weaker, more in- 
direct coupling, where either the matter dynamics in the 
disk lead to gravitational waves, or a global oscillation or 



FIG. 14: Same as Fig. [13] but we monitor the diagnostics at 
r = 3.95M (upper two panels) and 9.87M (lower two panels), 
respectively. 

resonance of whole spacetime exists, including both the 
exterior matter and the BH. 

In the first category where the wave emission is dom- 
inated by the motion of the fluid in Kerr spacetime, 
there are two possibilities. One is that inhomogeneities 
or radial motion in the disk exist, leading to the wave 
emission. These can be viewed, to a first approxima- 
tion, as small fluid fragments behaving approximately 
like test particles orbiting a Kerr BH. If we further ap- 
proximate the "particle" to be moving in a circular orbit 
then the characteristic frequency ^circular and the am- 
plitude /^circular from the quadrupole formula of gravita- 
tional waves can be written as 

^circular = TYl£l rQ , (3.12) 
rh+ circular = 7 X 10" 1 fl(M(j) 3/2 cos ut, (3.13) 

where Q ro is the orbital angular velocity of the particle 
at the radius of ISCO ro, and r is the distance from the 
source. The typical frequency for a circular orbit can be 
also applied for the scattering particle at the innermost 
radius ro- We show the relation between the orbital an- 
gular velocity at ISCO and the Kerr parameter in Fig.l23l 
In fact, when the newly formed BH is near extremal, the 
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FIG. 15: Spectra of m = 1 and m — 2 diagnostics throughout 
the collapse for model I. The top panel is monitored at the 
coordinate radius r = 1.02M, while the bottom panel r = 
2.05M. Note that the solid line and the dashed line denote 
m — 1 and m — 2. 



FIG. 16: Spectra of m = 1 and m — 2 diagnostics throughout 
the collapse for model II. The top panel is monitored at the 
coordinate radius r = 0.197M, while the bottom panel r = 
0.987M. Note that the solid line and the dashed line denote 
m — 1 and m — 2. 



typical frequency circular and the amplitude h+ circular 
of gravitational waves are 

MUJ circular M^ rQ « 0.5, 

* {Muj\ 3/2 

rh + circular ~ 2 X 10" 5 Af ( — J , 

assuming that the only contribution to the amplitude 
comes from the particle at ISCO (the ring mass is ~ 
3 x 10 -5 ), and that m = 1 is dominant (circular = ^r )- 
For ra = 1 we have that Mco MQ « 0.5. This idea 
is suggested from the spectra of gravitational waveform 
(Figs. n~UI and [20]). At every observed radius, there is a 
peak around Moo ^0.5, which might correspond to the 
motion of the particle at inner edge. However this idea 
cannot explain the enhancement of the quasi-periodic 
waves when the newly formed BH is almost extremal. 
We also note that the Fourier coefficients as shown in 
Figs. [TT]- HH show non-axisymmetries growing in the in- 
terior and moving out far too late (in retarded time) to 
cause the enhanced gravitational waves. 

Remaining with the idea that the gravitational waves 



are generated by the fluid motion, we note that the dy- 
namics of the matter as shown by Figs. [5] - [8] change 
after horizon formation. Initially all matter collapses in. 
However, once the quasi-stable disk forms in the interior 
there is still some infalling matter left, which "splashes" 
onto the disk. This is reminiscent of the centrifugal 
bounce seen in early core-collapse simulations with a sim- 
ple equation of state (e.g. the Type II multiple-bounce 
models of Ref. [3l|). As the rotation rate increases and 
the spacetime is near extreme Kerr, the chances of shocks 
appearing through some form of accretion are signifi- 
cantly increased, as shown as far back as Wilson [32| . 
After submission, we note that simulations of core col- 
lapse (|33[) have seen gravitational waves due to rapidly 
decelerated shocks. The shocks formed by the accreting 
matter could lead to wave-packet like bursts of gravita- 
tional waves. As these form in the outer layers of the disk 
and predominantly propagate out, their effects are never 
seen in the Fourier coefficients shown in Figs. [TT]- [TH It 
is possible to approximately correlate the shocks appear- 
ing in Figs. [2D and [22l with the start of the quasi-periodic 
waves in Figs. [TT] and [181 Rest-mass density snapshots 
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FIG. 17: Gravitational waveform measured with the Weyl 
scalar ^4 observed along the x-axis in the equatorial plane 
at r = 49.14M (top panel), r = 65.52M (middle panel) and 
98.29M (bottom panel) for model I. Note that the time at 
which the apparent horizon is first detected is t — 670M. 
Taking the wave propagation time from the source to the ob- 
server into account, the apparent horizon formation in the 
waveform is roughly just before the peak due to the burst. 
Initially a standard burst and ringdown signal are seen (com- 
pare [30j], especially their Fig. 18), but sustained gravitational 
wave emission indicate additional dynamics after BH forma- 
tion. 



in the meridional plane show that shock waves are gen- 
erated around t « 734M for model I and t « 828M for 
model II. Quasi-periodic waves in ^4 starts at t ~ 780M 
at r ^ 50M for model I, at t « 870M at r « 50M for 
model II. Since radiation propagates at speed of light, 
those times for each model roughly correspond with each 
other. However, it seems unlikely that these shocks could 
produce quasi-periodic waves of near constant amplitude 
as seen in model I, where the newly formed BH is not an 
extremal one (Fig. [T7|). 

In the second category, we are looking for a property 
of the full spacetime including both the BH and the mat- 
ter. As the Kerr BH contains the majority of the mass 
in the system, we would expect this to dominate. So ei- 
ther there must be some special property of quasi-normal 
modes of near extremal Kerr BHs, or there must be a cou- 
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FIG. 18: Gravitational waveform measured with the Weyl 
scalar ^4 observed along the x-axis in the equatorial plane 
at r = 47.36M (top panel), r = 63.15M (middle panel) and 
94.72M (bottom panel) for model II. Note that the time at 
which the apparent horizon is first detected is t = 770 M. 
Taking the wave propagation time from the source to the 
observer into account, the apparent horizon formation in the 
waveform is roughly just before the peak due to the burst. 



pling through the spacetime between the matter in the 
disk and the Kerr BH. Since we know that the imaginary 
part of the quasi-normal mode frequency goes to zero at 
extreme Kerr limit [34|], there are long-lived modes for 
near-extremal Kerr BHs, and hence small perturbations 
from the exterior matter may produce long-lived grav- 
itational wave signals such as those seen in Fig. [TT] for 
model I. It is not clear that this could easily lead to the 
strong amplification seen for model II in Fig. [T8l 

A final possibility that we consider in this category 
is a corotation resonance. In fact when the spin of the 
newly formed BH is very close to extreme Kerr, the am- 
plitude of the gravitational wave signal gradually grows 
after the quasinormal ringing. We also check the az- 
imuthal m modes of the rest mass density traced at the 
certain radius in the equatorial plane, and found that the 
m = 2 mode starts growing exponentially after the ring- 
down (Fig. [161 second top panel). One possible explana- 
tion for the exponential growth of the m = 2 mode at late 
times is the existence of corotation resonance of the newly 
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FIG. 19: Spectra of the real part of the Weyl scalar ^4 ob- 
served along the x-axis in the equatorial plane at r = 49.14M 
(top panel), r = 65.52M (middle panel) and 98.29M (bottom 
panel) for model I. 
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FIG. 20: Spectra of the Weyl scalar ^4 observed along the 
x-axis in the equatorial plane at r = 47.36M (top panel), 
r = 63.15M (middle panel) and 94.72M (bottom panel) for 
model II. 



formed disk triggered by the vibration of the hole. Note 
that the necessary condition to excite corotation reso- 
nance in fluid mechanics is cj c hr = mfl (e.g. Ref. [35[), 
where uo^ is the characteristic frequency of the system, 
Q the angular velocity. If the corotation resonance is 
triggered by quasinormal ringing, the necessary condition 
for triggering a corotation resonance is cj qnm = mft at a 
certain radius of the interior star. The real part of the 
quasi-normal mode frequency of I = m = 2 increases with 
increasing Kerr parameter, and becomes Mcc; qnm = 1 
for the extreme Kerr BEL The orbital velocity at IS CO 
also increases with increasing Kerr parameter, and be- 
comes Mft ro = (3 — v // 3)/2 = mcjparticie for the extreme 
Kerr BH (o; qnm < mo;p ar ticie)- Therefore, the condition 
w qnm = mcjparicie is satisfied for m = 2 for a certain Kerr 
parameter near the extremal limit. 



IV. CONCLUSION 

We have investigated the collapse of differentially ro- 
tating SMSs, especially focusing on the post BH forma- 
tion stage, by means of three dimensional hydrodynamic 



simulations in general relativity. We particularly focus 
on the onset of collapse to form a rapidly rotating BH as 
a final outcome. 

We have found that the qualitative results of the evo- 
lution for the mass and spin of the final BH and disk 
are quite similar to the estimates that can be computed 
from the equilibrium configuration when the estimated, 
final BH has ^(bh)/^(bh) ^ ^* This result suggests that 
in the absence of a nonaxisymmetric instability, the es- 
timate of the BH mass and the disk mass agree with a 
simple axisymmetric picture that the specific angular mo- 
mentum is conserved throughout the evolution, and the 
newly formed BH swallows the matter up to the radius 
of the ISCO. 

We have also found that a quasi-periodic wave oc- 
curs after the ringdown of a newly formed BH. As we 
would normally expect the ringdown waveform to damp, 
it seems likely that the cause of this waveform is due to 
the presence of the disk in some form. Furthermore, when 
the newly formed BH is sufficiently close to extreme Kerr 
with sufficient surrounding matter we have found that the 
wavesignal may be significantly amplified. 

We have discussed several possibilities for the origin 
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FIG. 21: Contours of the rest mass density in the meridional 
plane for model I. Snapshots are plotted at values of (£, p ma x) 
= i (713M, 1.87 x 10" 4 ), ii (734M, 1.24 x 10" 4 ), hi (755M, 
5.28 xlO -5 ) and iv (776M, 2.87 x 10" 5 ), with a cutoff density 
of pcut — 1.56 x 10 -13 , respectively. The contour lines denote 
rest mass densities p/p ma x = l0 _ ( 15_z ) d (i = 1, • • • , 14), where 
d = (log pmax — log p C ut)/15. Note that the apparent horizon 
exists after t = 671 M, and the coordinate radius of the ap- 
parent horizon in the equatorial plane are rh rz = i 0.295M, ii 
0.282M, hi 0.283M and iv 0.297M, respectively. Shocks can 
be found in the snapshots of I-ii. 
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FIG. 22: Contours of the rest mass density in the merid- 
ional plane for model II. Snapshots are plotted at values of 
(t, Pmax ) = i (808M, 3.15 x 1(T 4 ), ii (828M, 1.88 x 1(T 4 ), 
iii (849M, 1.05 x 10" 4 ) and iv (869M, 4.85 x 10" 5 ), with 
a cutoff density of p cu t = 1.56 x 10 -13 , respectively. The 
contour lines denote rest mass densities p/pmax = 10~( 15_ ^ d 
(i = 1, • • • ,14), where d = (logp max - log p cu t)/15. Note that 
the apparent horizon exists after t — 770M, and the coordi- 
nate radius of the apparent horizon in the equatorial plane 
are r hrz = i 0.173M, ii 0.167M, iii 0.170M and iv 0.180M, 
respectively. Shocks can be found in the snapshot of Il-ii. 



of these amplified waves. The most likely possibilities 
seem to be (a) corotation resonance between the disk 
and the BH, (b) long-lived gravitational waves from the 
near-extreme BH amplified by perturbations in the disk, 
or (c) shocks from the infalling, accreting matter. 

It may be possible to distinguish between some of 
these possibilities by further time integration, although 
at present the computational expense is prohibitive. 
The wave-packet or burst nature of waves generated by 
shocks, for example, should be possible to distinguish 



at later times and at larger radii of extraction. The slow 
power law decay (~ 1/t) in the extremal limit [36[ would, 
however, be difficult to accurately capture at present 
with this type of nonlinear numerical simulation. Al- 
ternatives to free evolution simulations such as the Fully 
Constrained Formalism of Refs. [3?], [3 should signifi- 
cantly improve the accuracy and hence may be capable 
of capturing such effects, but the specific gauge currently 
required by this formalism may make it difficult to get a 
numerically well-behaved simulation of the extreme cases 
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FIG. 23: Angular velocity of circular orbit at ISCO radius 
as a function of Kerr parameter. 



investigated here. 

Additional information about the horizon geometry 
could be found from the multipole structure of the dy- 
namical horizon (as defined in Ref. (39[; various numerical 
implementations are discussed in Refs. pol 14(1 l4ll l42| ) . 
In particular this may indicate whether a strong pertur- 
bation of the BH geometry is connected to the gravita- 
tional waves seen here. However, it seems likely, compar- 
ing in particular to the results of Refs. [4l|,[42[ where long 
term evolutions of "clean" binary black hole systems are 
considered, that better numerical accuracy and longer in- 
tegration times will be required in order to rule out any 
but the strongest of perturbations. 

Finally we discuss the detectability of gravitational 
waves from this system. The characteristic strength h 
and the frequency / of the burst and ringing are (e.g. 
Ref. 0) 



/burst 
^burst 

/qnm 

^QNM 



3 x 10~ 2 



io6M ^y /2 [ H z], 



M 



R ) 



2 x 10 



, , 1O 6 M , [Hz]; 



M 



6 x i -i9 ^w/M \ 1/a |< 2 x 10-^ [Hz] f 



/qnm 



1/2 



M 



V 10 6 M G 



1/2 m 



(4.1) 
(4.2) 
(4.3) 

(4.4) 



where d is the distance from the observer and A^gw 
is the total radiated energy. We set R/M = 1, a 
characteristic mean radius during BH formation. Since 
the main targets of LISA are gravitational radiation 
sources between 10 -4 and 10 _1 Hz, it is possible that 
LISA can search for the burst and quasi-normal ringing 
waves accompanying rotating SMS collapse and forma- 
tion of an SMBH. Moreover, the frequency of the quasi- 
periodic waves after the quasinormal ringing is quite sim- 
ilar to that of quasinormal ringing, so is in the detectable 
regime. Also the amplitude of the quasi-periodic wave is 
roughly 10% of the burst at the horizon, and hence this 
feature can also be seen in LISA. 

Although we focus on the case of the collapse of dif- 
ferentially rotating SMS to a SMBH in this paper, our 
model is also applicable to the collapse of population III 
stars which are also radiation pressure dominated. In this 
case, the mass range of the collapsing object becomes of 
the order of a few hundred solar masses. Therefore the 
most sensitive region for the detection of burst waves and 
quasinormal ringing waves becomes the order of hertz. 
Such gravitational waves might be seen in Deci-hertz in- 
terferometer gravitational wave observatory (DECIGO) 
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